module nem
    implicit none
    real(8), private, pointer, dimension(:) :: psi, psip, q
    real(8), private, pointer, dimension(:,:) :: phi
    real(8), private, pointer, dimension(:,:,:) :: outj
    real(8), private, pointer, dimension(:,:,:) :: inj
    real(8), private, pointer, dimension(:,:,:) :: coefm, coeff
    real(8), private, pointer, dimension(:,:,:) :: dfnem, dftemp
    integer, private :: dim
    integer, private, parameter :: fnem2 = 99, fnem4 = 199
    real(8), private :: jalpha (2)
   
contains
subroutine init_nem
    include "pspec.FI"
    include "xsec.FI"
    dim = totm_cmfd
    allocate(phi(dim,ng))
    allocate(outj(2,0:dim+1,ng))
    allocate(inj (2,0:dim+1,ng))
    allocate(dfnem(2,0:dim+1,ng))
    allocate(dftemp(2,0:dim+1,ng))
    allocate(coefm(4,dim,ng)) 
    allocate(coeff(4,dim,ng))   
    allocate(psi(dim))
    allocate(psip(dim))
    allocate(q(dim))
    open(fnem2, file = "RNEM2.txt", status = 'unknown')
    open(fnem4, file = "RNEM4.txt", status = 'unknown')
end subroutine

subroutine df_nem
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: l,i,m,g
    real(8) :: phm, pht, nphi
    real(8) :: A(2,2), x(2), r(2)
    real(8) :: ji(2), jo(2)
    dfnem = 1.
    dftemp = 1.

    do g=1, ng
        do l=1, dim
                nphi = hphi(l,g)
!            if (l.eq.1 .or. l.eq.dim) then
                m=2*l
!                ji(1) = hpj(m-2,g)
!                ji(2) = hpj(m+1,g)
!                jo(1) = jol(l,g,nphi,ji)
!                jo(2) = jor(l,g,nphi,ji)
!                do i=1, 2
!                    phm = 2*(ji(i)+jo(i))
!                    pht = 2*(hpj(m-1,g)+hpj(m,g))
!                    dfnem(i,l,g) = pht/phm
!                enddo
!            else
                A(1,1) = coefm(2,l,g)-1
                A(1,2) = coefm(3,l,g)
                A(2,1) = coefm(3,l,g)
                A(2,2) = coefm(2,l,g)-1
                r(1) = hj(2*l-1,g) - coefm(1,l,g)*nphi
                r(2) = hj(2*l  ,g) - coefm(1,l,g)*nphi
               
!                if (l .eq. 1) then
!                    A(:,1) = 0
!                    x(1) = 0
!                    !x(2) = r(1)/A(1,2)
!                    x(2) = r(2)/A(2,2)
!                elseif(l .eq. dim) then
!                    A(:,2) = 0
!                    x(2) = 0
!                    x(1) = r(1)/A(1,1)
!                    !x(1) = r(2)/A(2,1)
!                else
!                    call solve2x2(A,x,r)
!                endif
                call solve2x2(A,x,r)
                
                ji=x
                jo(1) = jol(l,g,nphi,ji)
                jo(2) = jor(l,g,nphi,ji)
                do i=1, 2
                    phm = 2*(ji(i)+jo(i))
                    pht = 2*(hpj(2*(i-1)+m-2,g)+hpj(2*(i-1)+m-1,g)) !hphis(2*(l-1)+i,g)
                    dftemp(i,l,g) = pht/phm
                enddo
                if (l.eq.1) jalpha(1) = ji(1)/jo(1)
                if (l.eq.dim) jalpha(2) = ji(2)/jo(2)
                inj(:,l,g) = ji
                outj(:,l,g) = jo
!            endif
        enddo
!        jcor(1) = hj(1,g)/(outj(1,1,g)-inj(1,1,g))
!        jcor(2) = hj(2*dim,g)/(outj(2,dim,g)-inj(2,dim,g))
    enddo
!            m=2*i-1
!            phm = 2*(hpj(m-1,g)+hpj(m,g))
!            pht = hphis(m,g)
!            dfnem(1,i,g) = pht/phm
!            m=2*i
!            phm = 2*(hpj(m,g)+hpj(m+1,g))
!            pht = hphis(m,g)
!            dfnem(2,i,g) = pht/phm

!            A(1,1) = coefm(2,l,g)-1
!            A(1,2) = coefm(3,l,g)
!            A(2,1) = coefm(3,l,g)
!            A(2,2) = coefm(2,l,g)-1
!            r(1) = hj(2*l-1,g) - coefm(1,l,g)*hphi(l,g)
!            r(2) = hj(2*l  ,g) - coefm(1,l,g)*hphi(l,g)
!            if (l .eq. 1) then
!                A(:,1) = 0
!                x(1) = 0
!                !x(2) = r(1)/A(1,2)
!                x(2) = r(2)/A(2,2)
!            elseif(l .eq. dim) then
!                A(:,2) = 0
!                x(2) = 0
!                x(1) = r(1)/A(1,1)
!                !x(1) = r(2)/A(2,1)
!            else
!                call solve2x2(A,x,r)
!            endif
!!            call solve2x2(A,x,r)
!            
!            ji=x
!            jo(1) = jol(l,g,hphi(l,g),ji)
!            jo(2) = jor(l,g,hphi(l,g),ji)
!            do i=1, 2
!                phm = 2*(ji(i)+jo(i))
!                pht = hphis(2*(l-1)+i,g)
!                dfnem(i,l,g) = pht/phm
!            enddo
!            inj(:,l,g) = ji
!            outj(:,l,g) = jo


!                if (l .eq. 1) then
!                    ji(1) = 0
!                    jo(1) = hj(1,g)
!                    ji(2) = (hj(2,g)-hj(1,g))/(coefm(2,l,g)-coefm(3,l,g)-1)
!                    nphi  = (hj(1,g)-coefm(3,l,g)*ji(2))/coefm(1,l,g)
!                    jo(2) = jor(l,g,nphi,ji)
!                    !x(2) = r(2)/A(2,2)
!                elseif(l .eq. dim) then
!!                    A(:,2) = 0
!!                    x(2) = 0
!!                    x(1) = r(1)/A(1,1)
!!                    !x(1) = r(2)/A(2,1)
!                    ji(2) = 0
!                    jo(2) = hj(2*l,g)
!                    ji(1) = (hj(2*l-1,g)-hj(2*l,g))/(coefm(2,l,g)-coefm(3,l,g)-1)
!                    nphi  = (hj(2*l,g)-coefm(3,l,g)*ji(1))/coefm(1,l,g)
!                    jo(1) = jol(l,g,nphi,ji)
!                    !x(2) = r(2)/A(2,2)
!                else
!                    ji(1) = 0.
!                    jo(1) = hj(1,g)
!                    A(1,1) = coefm(2,l,g)-1
!                    A(1,2) = coefm(3,l,g)
!                    A(2,1) = coefm(3,l,g)
!                    A(2,2) = coefm(2,l,g)-1
!                    r(1) = hj(2*l-1,g) - coefm(1,l,g)*nphi
!                    r(2) = hj(2*l  ,g) - coefm(1,l,g)*nphi
!                
!                    call solve2x2(A,x,r)
!                endif
!!                call solve2x2(A,x,r)
contains
subroutine solve2x2(A, x, r)
    real(8), intent(in) :: A(2,2), r(2)
    real(8), intent(out) :: x(2)
    real(8) :: detr
    detr = 1._8/(A(1,1)*A(2,2)-A(1,2)*A(2,1))
    x(1) = detr * (A(2,2)*r(1) - A(1,2)*r(2))
    x(2) = detr * (A(1,1)*r(2) - A(2,1)*r(1))
end subroutine

end subroutine

subroutine df_nem4
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: l,i,m,g
    real(8) :: phm, pht, nphi
    real(8) :: A(2,2), x(2), r(2)
    real(8) :: ji(2), jo(2)
    dfnem = 1.
    dftemp = 1.

    do g=1, ng
        do l=1, dim
                m=2*l
                call getpj(l,g,ji,jo)
                do i=1, 2
                    phm = 2*(ji(i)+jo(i))
                    dftemp(i,l,g) = 1._8/phm
                enddo
                if (l.eq.1) jalpha(1) = ji(1)/jo(1)
                if (l.eq.dim) jalpha(2) = ji(2)/jo(2)
                inj(:,l,g) = ji
                outj(:,l,g) = jo
        enddo
    enddo

contains
subroutine getpj(l,g,ji,jo)
    use LU_general
    use rng
    include "xsec.FI"
    include "hmgz.FI"
    integer, intent(in) :: l,g
    real(8), intent(out) :: ji(2), jo(2)
    real(8) :: A(4,4), x(4), r(4)
    real(8) :: c(4), M1, M2, N1, N2
    
    c = coefm(:,l,g)
    M1 = (hnuf(l,g)/hk+hmsct(g,g,l)-htot(l,g))*7/6
    M2 = htot(l,g)+42._8*hdif(l,g)/h_cmfd(l)**2-hmsct(g,g,l)-hnuf(l,g)/hk
    N1 = (-hnuf(l,g)/hk-hmsct(g,g,l)+htot(l,g))*3/2
    N2 = htot(l,g)+90._8*hdif(l,g)/h_cmfd(l)**2-hmsct(g,g,l)-hnuf(l,g)/hk
    A = 0
    A(1,1)=c(2)-1; A(2,1)=c(3);   A(3,1)=-c(4); A(4,1)=c(1)
    A(1,2)=c(3);   A(2,2)=c(2)-1; A(3,2)=c(4);  A(4,2)=c(1)
    A(1,3)=2;      A(2,3)=-2;     A(3,3)=M2/M1; A(4,3)=0
    A(1,4)=2;      A(2,4)=2;      A(3,4)=0;     A(4,4)=N2/N1
    r(1) = hj(2*l-1,g)-c(1)*hphi(l,g); r(2) = hj(2*l,g)-c(1)*hphi(l,g)
    r(3) = hj(2*l,g)-hj(2*l-1,g);      r(4) = hphi(l,g)-(hj(2*l,g)+hj(2*l-1,g))
    call LU_factorize(4,A)
    call solve_LU(4,A,r,x)
    ji = x(1:2)
    jo(1) = ji(1)+hj(2*l-1,g)
    jo(2) = ji(2)+hj(2*l,g)
end subroutine

subroutine getpj_old(l,g,ji,jo)
    use LU_general
    use rng
    include "xsec.FI"
    include "hmgz.FI"
    integer, intent(in) :: l,g
    real(8), intent(out) :: ji(2), jo(2)
    real(8) :: A(5,5), x(5), r(5)
    real(8) :: c(4), xi, zeta
    integer :: i, j
    
    c = coefm(:,l,g)
    A = 0
    A(1,1)=c(2); A(2,1)=-1;   A(3,1)=c(3); A(5,1)=-c(4)
    A(1,2)=c(3); A(3,2)=c(2); A(4,2)=-1;   A(5,2)=c(4)
    A(1,3)=-1;   A(2,3)=1
    A(3,4)=-1;   A(4,4)=1
    xi = (hnuf(l,g)/hk+hmsct(g,g,l)-htot(l,g))*7/6
    zeta = htot(l,g)+42._8*hdif(l,g)/h_cmfd(l)**2-hmsct(g,g,l)-hnuf(l,g)/hk
    A(1,5)=xi; A(2,5)=xi; A(3,5)=-xi; A(4,5)=-xi; A(5,5)=zeta
    r(1) = -c(1)*hphi(l,g); r(2) = -c(1)*hphi(l,g)
    r(3) = hj(2*l-1,g); r(4) = hj(2*l,g); r(5) = 0
    call LU_factorize(5,A)
    call solve_LU(5,A,r,x)
    ji(1) = x(1); ji(2) = x(3)
    jo(1) = x(2); jo(2) = x(4)
end subroutine

end subroutine

function jor(l, g, phi, ji)
    integer, intent(in) :: l, g
    real(8), intent(in) :: phi, ji(2)
    real(8) :: jor
    jor = coefm(1,l,g)*phi + coefm(3,l,g)*ji(1) + coefm(2,l,g)*ji(2)
end function

function jol(l, g, phi, ji)
    integer, intent(in) :: l, g
    real(8), intent(in) :: phi, ji(2)
    real(8) :: jol
    jol = coefm(1,l,g)*phi + coefm(2,l,g)*ji(1) + coefm(3,l,g)*ji(2)
end function

subroutine solve_nem_df
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    real(8) :: k, err
    real(8) :: jil, jir, jor, jol, leak, alpha
    integer :: i,l,g,iter
    logical, save :: run = .false.
    call make_coefm
    call df_nem
!    if (run .eq. .false.) then
!        dftemp = dfnem
!        run = .true.
!    endif
    dfnem = dftemp
    phi = hphi
    k=hk

!    inj = 0.25
!    outj = 0.25
    outj(:,0,:) = 0
    outj(:,dim+1,:) = 0
    inj(:,0,:) = 0
    inj(:,dim+1,:) = 0
!    inj(1,1,:) = 0.
!    inj(2,dim,:) = 0.
    
!    dfnem = dftemp
    iter = 0
    err = 1.
    call makeFS
    do while (err.gt.1e-6)  
        iter = iter+1
        do g=1, ng
            call make_source(g,k)
            
            ! partial current update
            do i=1, dim
                l = redblack(i,dim)
                ! outgoing partial current update
                outj(1,l,g) = coefm(1,l,g)*phi(l,g) + coefm(2,l,g)*inj(1,l,g) + coefm(3,l,g)*inj(2,l,g)
                outj(2,l,g) = coefm(1,l,g)*phi(l,g) + coefm(3,l,g)*inj(1,l,g) + coefm(2,l,g)*inj(2,l,g)
                ! incoming partial current update
                alpha = 0.5_8*(1-dfnem(1,l,g)/dfnem(2,l-1,g))
                inj(2,l-1,g) = (1-alpha)*outj(1,l,g)-alpha*inj(1,l,g)
                alpha = 0.5_8*(1-dfnem(2,l,g)/dfnem(1,l+1,g))
                inj(1,l+1,g) = (1-alpha)*outj(2,l,g)-alpha*inj(2,l,g)
            end do
!            jalpha(1) = 1-hj(1,g)/outj(1,1,g)
!            jalpha(2) = 1-hj(2*dim,g)/outj(2,dim,g)
!            jalpha = 0.5_8
            inj(1,1,g) = jalpha(1) * outj(1,1,g)
            inj(2,dim,g) = jalpha(2) * outj(2,dim,g)
            
!            inj(1,1,g) = outj(1,1,g)-hj(1,g)            
!            inj(2,dim,g) = outj(2,dim,g)-hj(2*dim,g)
            ! node average flux update
            do i=1, dim
                l = redblack(i,dim)
                leak = (outj(2,l,g)-inj(2,l,g)) - (inj(1,l,g)-outj(1,l,g))
                phi(l,g) = (q(l)-leak/h_cmfd(l))/htot(l,g)
            enddo
        end do
        call estimateK(k, err)
        if (iter==30000) exit
    end do
!    call init_nem
    call makeFS
    phi = phi/sum(psi*h_cmfd(1))
    write(fnem2, '(2000(e11.5,1x))') k,(phi(i,1),i=1,dim)
end subroutine

subroutine solve_nem_df4
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    real(8) :: k, err
    real(8) :: jil, jir, jor, jol, leak, alpha
    integer :: i,l,g,iter
    logical, save :: run = .false.
    call make_coefm
    call df_nem4
!    if (run .eq. .false.) then
!        dftemp = dfnem
!        run = .true.
!    endif
    dfnem = dftemp
    phi = hphi
    k=hk

!    inj = 0.25
!    outj = 0.25
    outj(:,0,:) = 0
    outj(:,dim+1,:) = 0
    inj(:,0,:) = 0
    inj(:,dim+1,:) = 0
!    inj(1,1,:) = 0.
!    inj(2,dim,:) = 0.
    
!    dfnem = dftemp
    iter = 0
    err = 1.
    call makeFS
    do while (err.gt.1e-6)
        iter = iter+1
        do g=1, ng
            call make_source(g,k)
            
            ! partial current update
            do i=1, dim
                l = redblack(i,dim)
                ! outgoing partial current update
                coeff(1,l,g) = coef1(l,g)
                coeff(2,l,g) = coef2(l,g)
                coeff(3,l,g) = coef3(k,l,g)
                coeff(4,l,g) = coef4(k,l,g)
                outj(1,l,g) = coefm(1,l,g)*(phi(l,g)+coeff(4,l,g)) + coefm(2,l,g)*inj(1,l,g) + coefm(3,l,g)*inj(2,l,g) - coefm(4,l,g)*coeff(3,l,g)
                outj(2,l,g) = coefm(1,l,g)*(phi(l,g)+coeff(4,l,g)) + coefm(3,l,g)*inj(1,l,g) + coefm(2,l,g)*inj(2,l,g) + coefm(4,l,g)*coeff(3,l,g)
                ! incoming partial current update
                alpha = 0.5_8*(1-dfnem(1,l,g)/dfnem(2,l-1,g))
                inj(2,l-1,g) = (1-alpha)*outj(1,l,g)-alpha*inj(1,l,g)
                alpha = 0.5_8*(1-dfnem(2,l,g)/dfnem(1,l+1,g))
                inj(1,l+1,g) = (1-alpha)*outj(2,l,g)-alpha*inj(2,l,g)
            end do
            inj(1,1,g) = jalpha(1) * outj(1,1,g)
            inj(2,dim,g) = jalpha(2) * outj(2,dim,g)
            
            ! node average flux update
            do i=1, dim
                l = redblack(i,dim)
                leak = (outj(2,l,g)-inj(2,l,g)) - (inj(1,l,g)-outj(1,l,g))
                phi(l,g) = (q(l)-leak/h_cmfd(l))/htot(l,g)
            enddo
        end do
        call estimateK(k, err)
        if (iter==30000) exit
    end do
!    call init_nem
    call makeFS
    phi = phi/sum(psi*h_cmfd(1))
    write(fnem4, '(2000(e11.5,1x))') k,(phi(i,1),i=1,dim)
end subroutine

subroutine solve_nem4
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    real(8) :: k, err
    real(8) :: jil, jir, jor, jol, leak
    integer :: i,l,g,iter
    phi = 1._8
    outj = 1._8/4
    outj(:,0,:) = 0
    outj(:,dim+1,:) = 0
    inj = outj
    inj(1,1,:) =0
    inj(2,dim,:) = 0
    k=1
    iter = 0
    err = 1.
    
    call make_coefm
    
    call makeFS
    do while (err.gt.1e-10)
        iter = iter+1
        do g=1, ng
            call make_source(g,k)
            
            ! partial current update
            do i=1, dim
                l = redblack(i,dim)
                coeff(1,l,g) = coef1(l,g)
                coeff(2,l,g) = coef2(l,g)
                coeff(3,l,g) = coef3(k,l,g)
                coeff(4,l,g) = coef4(k,l,g)
                outj(1,l,g) = coefm(1,l,g)*(phi(l,g)+coeff(4,l,g)) + coefm(2,l,g)*inj(1,l,g) + coefm(3,l,g)*inj(2,l,g) - coefm(4,l,g)*coeff(3,l,g)
                outj(2,l,g) = coefm(1,l,g)*(phi(l,g)+coeff(4,l,g)) + coefm(3,l,g)*inj(1,l,g) + coefm(2,l,g)*inj(2,l,g) + coefm(4,l,g)*coeff(3,l,g)
                inj(2,l-1,g) = outj(1,l,g)
                inj(1,l+1,g) = outj(2,l,g)
            end do
            
            ! node average flux update
            do i=1, dim
                l = redblack(i,dim)
                jil = outj(2,l-1,g)
                jir = outj(1,l+1,g)
                jor = outj(2,l,g)
                jol = outj(1,l,g)
                leak = (jor-jir) - (jil-jol)
                phi(l,g) = (q(l)-leak/h_cmfd(l))/htot(l,g)
            enddo
        end do
        call estimateK(k, err)
        if (iter==30000) exit
    end do
end subroutine

function coef1(l,g)
    integer, intent(in) :: l, g
    real(8) :: coef1
    coef1 = (inj(2,l,g)+outj(2,l,g))-(inj(1,l,g)+outj(1,l,g))
end function
function coef2(l,g)
    integer, intent(in) :: l, g
    real(8) :: coef2
    coef2 = phi(l,g)-(inj(2,l,g)+outj(2,l,g)+inj(1,l,g)+outj(1,l,g))
end function

function coef3(k,l,g)
    include "hmgz.FI"
    include "xsec.FI"
    include "pspec.FI"
    real(8), intent(in) :: k
    integer, intent(in) :: l,g
    real(8) :: coef3
    real(8) :: rhs, lhs
    integer :: gp
    rhs = -htot(l,g)*coeff(1,l,g)
    do gp=1, ng
        rhs = rhs+(xchi(1,g)/k*hnuf(l,gp)*coeff(1,l,gp)+hmSct(g,gp,l)*coeff(1,l,gp))
    enddo
    rhs = rhs*7/6
    
    lhs = htot(l,g)+hdif(l,g)/h_cmfd(l)**2*42
    do gp=1, ng
        lhs = lhs - xchi(1,g)/k*hnuf(l,gp)
    enddo    
    lhs = lhs - hmSct(g,g,l)
    coef3 = rhs/lhs
end function

function coef4(k,l,g)
    include "hmgz.FI"
    include "xsec.FI"
    include "pspec.FI"
    real(8), intent(in) :: k
    integer, intent(in) :: l,g
    real(8) :: coef4
    real(8) :: rhs, lhs
    integer :: gp
    rhs = htot(l,g)*coeff(2,l,g)
    do gp=1, ng
        rhs = rhs - (xchi(1,g)/k*hnuf(l,gp)*coeff(2,l,gp)+hmSct(g,gp,l)*coeff(2,l,gp))
    enddo
    rhs = rhs*3/2
    
    lhs = htot(l,g)+hdif(l,g)/h_cmfd(l)**2*90
    do gp=1, ng
        lhs = lhs - xchi(1,g)/k*hnuf(l,gp)
    enddo    
    lhs = lhs - hmSct(g,g,l)
    coef4 = rhs/lhs 
end function

subroutine solve_nem
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    real(8) :: k, err
    real(8) :: jil, jir, jor, jol, leak
    integer :: i,l,g,iter
    phi = 1._8
    outj = 1._8/4
    outj(:,0,:) = 0
    outj(:,dim+1,:) = 0
    inj = outj
    inj(1,1,:) =0
    inj(2,dim,:) = 0
    k=1
    iter = 0
    err = 1.
    
    call make_coefm
    
    call makeFS
    do while (err.gt.1e-10)
        iter = iter+1
        do g=1, ng
            call make_source(g,k)
            
            ! partial current update
            do i=1, dim
                l = redblack(i,dim)
                outj(1,l,g) = coefm(1,l,g)*phi(l,g) + coefm(2,l,g)*inj(1,l,g) + coefm(3,l,g)*inj(2,l,g)
                outj(2,l,g) = coefm(1,l,g)*phi(l,g) + coefm(3,l,g)*inj(1,l,g) + coefm(2,l,g)*inj(2,l,g)
                inj(2,l-1,g) = outj(1,l,g)
                inj(1,l+1,g) = outj(2,l,g)
            end do
            
            ! node average flux update
            do i=1, dim
                l = redblack(i,dim)
                jil = outj(2,l-1,g)
                jir = outj(1,l+1,g)
                jor = outj(2,l,g)
                jol = outj(1,l,g)
                leak = (jor-jir) - (jil-jol)
                phi(l,g) = (q(l)-leak/h_cmfd(l))/htot(l,g)
            enddo
        end do
        call estimateK(k, err)
        if (iter==30000) exit
    end do
end subroutine



subroutine estimateK(k, err)
    real(8), intent(out) :: k, err
    real(8) :: upper, lower
    integer :: i
    ! evaluate k    
    call makeFS                                        
    upper = 0.
    lower = 0.
    do i=1, dim                                                
        upper = upper + psi(i)*psi(i)
        lower = lower + psi(i)*psip(i)
    end do
    k = k*upper/lower
    
    ! evaluate err 
    err = 0.
    lower = sqrt(upper)
    upper = 0
    do i=1, dim
        upper = upper + (psi(i)- psip(i))*(psi(i)- psip(i))
    end do
    upper = sqrt(upper)
    err = upper/lower
end subroutine

function redblack(i, dim) result(l)
    integer, intent(in) :: i, dim
    integer, save :: base = -1
    integer, save :: cnt = 0  
    integer :: l
    l = base + 2*i
    if (l>dim) then
        base = -(dim+mod(dim,2))
        l = base + 2*i
    endif
    cnt = cnt + 1
    if (cnt .eq. dim) then 
        cnt = 0
        base = -1
    endif
end function

    
subroutine make_coefm
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i,g
    real(8) :: beta, beta2
    do g=1, ng
        do i=1, dim
            beta = hdif(i,g)/h_cmfd(i)
            beta2 = beta*beta
            coefm(1,i,g) = 6._8*beta/(1 + 12*beta)
            coefm(2,i,g) = (1._8 - 48*beta2)/(1 + 16*beta + 48*beta2)
            coefm(3,i,g) = (-8._8*beta)/(1 + 16*beta + 48*beta2)
            coefm(4,i,g) = 6._8*beta/(1 + 4*beta)
        end do
    end do
end subroutine

subroutine make_source(g,k)
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer, intent(in) :: g
    real(8), intent(in) :: k
    integer :: i, gp
    do i=1, dim
        do gp=1, ng
            q(i) = psi(i)/k+hmSct(g,gp,i)*phi(i,gp)
        end do
    end do
end subroutine

subroutine makeFS
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i,g
    psip = psi
    psi = 0.
    do g=1, ng
        do i=1, dim
            psi(i) = psi(i) + hNuF(i,g) * phi(i, g)
        end do
    end do
end subroutine
end module
